\documentclass[a4paper]{article}
\usepackage[english]{babel}
\usepackage[T1]{fontenc}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amsthm}
\usepackage{graphicx}
\usepackage{algorithmic}
\usepackage{algorithm}
\usepackage{amssymb}
\usepackage{url}

\renewcommand{\algorithmicrequire}{\textbf{Input:}}
\renewcommand{\algorithmicensure}{\textbf{Output:}}
\usepackage[table]{xcolor}
\usepackage{color}
\usepackage{listings}
\usepackage{courier}
\lstset{ %
language=Java,                % choose the language of the code
basicstyle=\tiny,       % the size of the fonts that are used for the code
numbers=left,                   % where to put the line-numbers
numberstyle=\tiny,      % the size of the fonts that are used for the line-numbers
stepnumber=1,                   % the step between two line-numbers. If it is 1 each line will be numbered
numbersep=5pt,                  % how far the line-numbers are from the code
backgroundcolor=\color{white},  % choose the background color. You must add \usepackage{color}
showspaces=false,               % show spaces adding particular underscores
showstringspaces=false,         % underline spaces within strings
showtabs=false,                 % show tabs within strings adding particular underscores
frame=single,   		% adds a frame around the code
tabsize=2,  		% sets default tabsize to 2 spaces
captionpos=b,   		% sets the caption-position to bottom
breaklines=true,    	% sets automatic line breaking
breakatwhitespace=false,    % sets if automatic breaks should only happen at whitespace
escapeinside={\%}{)}          % if you want to add a comment within your code
}



\title{Advanced Algorithms \\ Project 1}
\author{David Falk, Klaus Nicosia, Marcus Öberg\\
dfa@kth.se, nicosia@kth.se, marcusob@kth.se}












\begin{document}
\maketitle
\thispagestyle{empty}
\newpage
\begin{abstract}
In this project we discuss the problem of finding the prime factors of integers. We then describe a couple of different methods of solving this problem and gives a presentation of how we chose to implement some of these in java and C++. In particular we describes the \textsc{Pollard-Rho} and the \textsc{Pollard-Brent} algorithms. In the results we saw that the \textsc{Pollard-Brent}  proved to be faster then \textsc{Pollard-Rho} algorithm. We come to the conclusion that the results we got is mainly due to the fact that \textsc{Pollard-Brent} calls GCD far less often then the \textsc{Pollard-Rho} algorithm.
\end{abstract}
\thispagestyle{empty}
\newpage
\tableofcontents
\thispagestyle{empty}
\newpage
\setcounter{page}{4}




\section{Introduction}
In this project we were supposed to find solutions to the factorization problem. In other words implement an algorithm that could find all prime factors of big integers as fast as possible. We have studied a couple of different algorithm and chosen to implement three different kinds of  Pollard algorithms. We have then evaluated the result with our own tests and through the validation system \textsc{Kattis}. In the end we discuss the result and how it could have been improved.
\section{Problem}

 Given an integer $N$ return all the prime factors of $N$.
\section{Possible solutions}

\subsection{General Principle}
All of these possible solutions only calculate one factor (not necessarily prime). Let's call this factor $f_1$. The other factor is calculated by $f_2 = N/f_1$. Then the factors are calculated recursively for $f_1$ and $f_2$. The recursion ends when a number is prime. The number $x$ is judged to be prime by using a function called \texttt{isPrime(x)}. If this function returns \texttt{true} then the number is a prime and no factorization is done.

If $N$ is even all factors of two are removed before using the factorization algorithm.

\subsubsection*{Psuedocode}
\begin{algorithm}
\caption{\textsc{Factor-Integer}(n)}
\label{Factor-Integer}
\begin{algorithmic}
\REQUIRE Number $n$ to factor
\ENSURE List with all factors
\STATE $list \Leftarrow \emptyset$
\IF{isPrime($n$)}
	\STATE $list.add(n)$
	\RETURN $list$
\ENDIF
\IF{$n$ is even}
	\STATE $list.add(2)$
	\STATE $newlist \Leftarrow$ \textsc{Factor-Integer}($n/2$)
	\STATE $list.add(newlist)$
\ELSE
	\STATE $p \Leftarrow $ \textsc{Factorization-Algorithm}(n)
	\STATE $newlist \Leftarrow$ \textsc{Factor-Integer}($p$)
	\STATE $list.add(newlist)$
	\STATE $newlist \Leftarrow$ \textsc{Factor-Integer}($n/p$)
	\STATE $list.add(newlist)$
\ENDIF
\RETURN $list$
\end{algorithmic}
\end{algorithm}

\subsection{Na\"{\i}ve Solution}
The na\"{\i}ve solution to solve the factorization problem is to start with $p=2$ and test if $p | N$ for all primes $p\leq \sqrt{N}$. If the prime divides $N$ ($p | N$) it is a prime factor of $N$. This solution is obviously very slow and (in worst case) all primes between $2$ and $\sqrt{N}$ have to be calculated.
 
\subsection{Fermat's Factorization Method}
 Fermat's factorization method factorizes $N$ by defining $a$ and $b$
 \begin{equation}
 \label{eq:fermat1}
 N = a^2 - b^2
 \end{equation}
 This is the same as
 \begin{equation}
 \label{eq:fermat2}
 N = (a+b)(a-b)
 \end{equation}

Every positive odd non-prime integer can be stated as in \eqref{eq:fermat1} by writing it as $N = xy$ and
 \begin{equation}
 \begin{array}{lcl} x = (a+b) \\ y = (a-b) \end{array}
 \end{equation}

First we create an integer $q = \lceil\sqrt{N}\rceil$.
Calculate $c = q^2-N$ and if $c$ is not a square number then $c$ is recalculated using $q = q+1$. When $c$ is a square number we have
 \begin{equation}
 \begin{array}{lcl} q = a \\ c = b^2 \end{array}
 \end{equation}
 which means that
 \begin{equation}
 \begin{array}{lcl} y = q-\sqrt{c} \\ x = N/y \end{array}
 \end{equation}

When we have found $x$ and $y$ we can now recursively use Fermat's factorization to calculate their factors.

\subsubsection*{Psuedocode}
\begin{algorithm}
\caption{\textsc{Fermats-Factorization-Algorithm}(n)}
\label{Fermat's-Factorization-Algorithm}
\begin{algorithmic}
\REQUIRE Odd number $n$ to factor
\ENSURE Factor $f$ of $n$
\STATE $q \Leftarrow$\textsc{CEIL}(\textsc{SQRT}$(n)$)
\STATE $c \Leftarrow q^2 - n$
\WHILE{$c$ is not a square}
	\STATE $q \Leftarrow q + 1$
	\STATE $c \Leftarrow q^2 - n$
\ENDWHILE
\RETURN $q-$\textsc{SQRT}$(c)$
\end{algorithmic}
\end{algorithm}

\subsection{Pollard-Rho}
The core idea in the in the Pollard-Rho algorithm is based on observations from the birthday problem (see section \ref{birthday}) and Floyd's cycle-finding algorithm (see section \ref{floyd}). From the birthday problem is can be calculated that two numbers x and y are congruent modulo p with a probability of 0,5 after $1,77\sqrt{p}$ random numbers have been chosen.  If n is the number we are trying to factor and if p is a factor of n, then 
\begin{equation}
p \leq gcd(x-y,n) \leq n
\end{equation}
since p divides both $x-y$ and $n$.

With the help of a function on the form $f(x) = kx^y+c\; mod\;n$ you can then exploit Floyd's cycle-finding algorithm (because of mod n, f(x) is in the ring $\mathbb{Z}_n$ and thus we can use Floyd's cycle-finding algorithm) to identify the numbers x and y.

\subsubsection{The Birthday Problem}\label{birthday}
The birthday problem or birthday paradox concerns the probability that two persons in a group with size n has the same birthday. From the pigeonhole principle, this states that if you put k items into m pigeonholes where k > m at least one hole will contain more then one item, you can conclude that with 366 persons (based on that there is 365 days in a year) there is 100\% probability that two persons will have the same birthday. For the probability calculation we assume that all persons are born the same year and that this year has$k= 365$ days, further we assume that all birthdays are independent, meaning that there are no twins or seasonal dependencies.

 The probability that one person is born at any day is $365/365$ or 100\%. So for two people to not have the same birthday the probability would be $365/365 \times 364/365$ (according to the multiplication principle). But if you consider the fact that we are looking for any two persons in an n sized group that share birthday we can conclude that if n is three the probability would be $365/365 \times 364/365 \times 363/364$ since the third person would not be allowed to share birthday with either the first person or the second person. If you consider the fact that $\bar{P}(A) = 1 - P(A)$ and that the probability is one when n > 365 (pigeonhole principle), you would get the following probability that two persons doesn't share birthday in an n sized group.

 
\begin{equation}
 \begin{array}{lcl}
\bar{P}(A) = 1 \times (1 - \frac{1}{365}) \times (1 - \frac{2}{365}) \times \dots\times (1 -\frac{n-1}{365})
\\
\\

= \frac{365 \times 364 \times \dots \times (365 - n + 1)}{365 ^n}
\\
\\
= \frac{365!}{365^n(365 - n)!} = \frac{n! \cdot \binom{365}{n} }{365^n}
\end{array}
\end{equation}

From this you could calculate that if n is 23 there would be about 50\% chance that two persons doesn't share birthday (and also that two persons share birthday) or that if n is 57 there is about 1\% chance that two persons doesn't share birthday (and 99\% that two persons do share birthday). The Pollard-Rho algorithm uses the fact that there is a very small probability loss when you decrease the size of n when n is close to k.


\subsubsection{Floyd's Cycle-Finding Algorithm}\label{floyd}
Floyd's cycle-finding algorithm is used to identify a cycle in given function call with input that maps to itself. For example given a set $S=\{1,2,3,4\}$ and a mapping of f(x)


\begin{equation}
\begin{array}{lcl}
f(1) = 2\\
f(2) = 3\\
f(3) = 4\\
f(4) = 3\\ 
\end{array}
\end{equation}

If you start at $x=1$ you would get the sequence $1, 2, 3, 4, 3, 4,\dots$ After you have reached $x=3$ you are stuck in a cycle. Floyd's way of solving this problem is sometimes referred to as the "tortoise and the hare algorithm". The analogy with the tortoise and the hare is actually quite useful and goes like this. A hare and a tortoise races through a course similar to the one in the example above. If the hare moves at double the speed  of the tortoise and they move in a discrete manner (in steps as opposed to linear) it is guaranteed that they will be at the same spot at some point in time. The reason that this is true is obvious since they move in a cycle.

For any integer $i > \mu$ where $ \mu$ is the first integer in the cycle and $k \geq 0$, $x_i = x_{i+k \lambda}$ where $\lambda$ is the length of the cycle to be found. This stated the obvious fact that if you are at any point in the cycle and move one or more full laps in the cycle you will end up in the same spot.From this it also follow that whenever $i = k\lambda \geq \mu$, $x_i = x_{2i}$ . This means that it is sufficient for the algorithm to check if the values of the two pointer (hare and tortoise) to be equal when the second pointer has moved two times as long as the first. This is the reason that the hare moves double the speed of the tortoise.



\subsubsection*{Psuedocode}
\begin{algorithm}
\caption{\textsc{Pollard-Rho}(n)}
\label{Pollard-Rho}
\begin{algorithmic}
\REQUIRE Odd number $n$ to factor
\ENSURE Factor $f$ of $n$
\STATE $x \Leftarrow$ Random integer in the range $[0,n-1]$
\STATE $y \Leftarrow x$
\STATE $d \Leftarrow 1$
\WHILE{$d = 1$}
	\STATE $x \Leftarrow f(x)$
	\STATE $y \Leftarrow f(f(y))$
	\STATE $d \Leftarrow $\textsc{GCD}$(n,\textsc{ABS}(x-y))$
\ENDWHILE
\RETURN $d$
\end{algorithmic}
\end{algorithm}

\subsection{Pollard-Rho with Brent's Cycle Detection}
Richard Brent improved \textsc{Pollard-Rho} with a better cycle detection. This algorithm uses the Brent's cycle-finding algorithm \ref{brent-cycle} opposed to the Floyd's cycle-finding algorithm \ref{floyd}.

\subsubsection{Brent's Cycle-Finding Algorithm}
\label{brent-cycle}
Brent's cycle-finding algorithm is similar to Floyd's cycle-finding algorithm \ref{floyd} in that it uses two values (tortoise and hare). Lets call these values $x$ and $y$ instead where $x_i$ is the $i$:th value of $x$.
For example given a set $S=\{1,2,3,4,5,6,7,8,9\}$
\begin{equation}
\begin{array}{lcl}
f(1) = 2\\
f(2) = 3\\
f(3) = 4\\
f(4) = 5\\
f(5) = 6\\
f(6) = 7\\
f(7) = 8\\
f(8) = 7\\
f(9) = 8\\
\end{array}
\end{equation}
If you start at $x=1$ you would get the sequence $1, 2, 3, 4, 5, 6, 7, 8, 7 ...$. After you have reached $x=7$ you are stuck in a cycle. Brent's way of finding the cycle is to set $y_i=x_{2^i}$ for all $i$. The reasons this works is the same as in Floyd's cycle finding algorithm with the only difference being how often one of the values are updated. Since Brent uses $x$ to calculate $y$ the function only needs to be called once in each step instead of three times in Floyd's.

\subsubsection*{Psuedocode}
\begin{algorithm}
\caption{\textsc{Pollard-Rho-Brent}(n)}
\label{Pollard-Rho-Brent}
\begin{algorithmic}
\REQUIRE Odd number $n$ to factor
\ENSURE Factor $f$ of $n$
\STATE $i \Leftarrow 1$
\STATE $k \Leftarrow 2$
\STATE $x \Leftarrow$ Random integer in the range $[0,n-1]$
\STATE $y \Leftarrow x$
\STATE $d \Leftarrow 1$
\WHILE{$true$}
	\STATE $i \Leftarrow i+1$
	\STATE $x \Leftarrow f(x)$
	\STATE $d \Leftarrow$\textsc{GCD}($n$,\textsc{ABS}$(x-y)$)
	\IF{$d = abs(y-x)$}
		\RETURN d
	\ENDIF
	\IF{$i = k$}
		\STATE $y \Leftarrow x$
		\STATE $k \Leftarrow 2k$
	\ENDIF
\ENDWHILE
\RETURN $d$
\end{algorithmic}
\end{algorithm}

\subsection{Pollard-Brent}
The \textsc{GCD} call is the slowest calculation in the different Pollard algorithms above. Brent used the fact that if \textsc{GCD}$(a, n) > 1$, then \textsc{GCD}$(ab, n) > 1$, were a and b are integers \cite{wikipedia_pollard}. In particular if \textsc{GCD}$(|x-y|, n) > 1$, then \textsc{GCD}$(|x-y|b, n) > 1$, so that instead of calculating \textsc{GCD}$(|x-y|, n)$ in every step, we instead calculate one \textsc{GCD}$(q, n)$, were $q$ is $m$ consecutive products of $|x-y|$ modulo $n$ \cite{brent}. So most of the \textsc{GCD} calculations are reduced to a multiplication modulo $n$.

If $m$ increases the probability to find a factor decreases, so $m$ should not be too large. Should this part of the algorithm fail we should backtrack to the $y$ value before \textsc{GCD}$(q, n)$ has been computed, to minimize this risk. This value is stored in $ys$ and can be used in Pollard's original algorithm.

From the sequence $y_1, y_2, y_3, \dots$ a corresponding sequence of $x_1, x_2, x_3, \dots$ is generated. For $x_{i} = y_{j}$, where $i$ is a postive integer, the index $j$ follows a similar, but more complex, series compared to that of Brent's cycle detection.

 \begin{algorithm}
\caption{\textsc{Pollard-Brent}(n)}
\label{Pollard-Brent}
\begin{algorithmic}
\REQUIRE Odd number $n$ to factor
\ENSURE Factor $f$ of $n$, $-1$ if failed
\STATE $r \Leftarrow 1$
\STATE $m \Leftarrow bitlength(n)$
\STATE $y \Leftarrow$\textsc{Random}(0, $n-1$)
\STATE $q \Leftarrow 1$
\STATE $g \Leftarrow 1$
\STATE $ys \Leftarrow y$
\WHILE{$g = 1$}
	\STATE $x \Leftarrow y$
	\FOR{$i = 0$ \TO $i = r-1$}
		\STATE $y \Leftarrow f(y)$
	\ENDFOR
	\STATE $k \Leftarrow 0$
	\WHILE{$k<r$ \AND $g = 1$}
		\STATE $ys \Leftarrow y$
		\FOR{$i = 0$ \TO $i = min(m$, $r-k)-1$}
			\STATE $y \Leftarrow f(y)$
			\STATE $q \Leftarrow q(x-y) \; mod \; n$
		\ENDFOR
		\STATE $g \Leftarrow$ \textsc{GCD}($q$, $n$)
		\STATE $k \Leftarrow k+m$
	\ENDWHILE
	\STATE $r \Leftarrow 2r$
\ENDWHILE
\IF{$g=n$}
	\STATE $ys \Leftarrow f(ys)$
	\STATE $g \Leftarrow gcd(n,x-ys)$
	\WHILE{$g=1$}
		\STATE $ys \Leftarrow f(ys)$
		\STATE $g \Leftarrow gcd(n,x-ys)$
	\ENDWHILE
\ENDIF
\IF{$g=n$}
	\RETURN $-1$
\ENDIF
\RETURN $g$
\end{algorithmic}
\end{algorithm}
 
\newpage 
 
 \subsection{Quadratic Sieve}

The quadratic sieve is the fastest known algorithm for big integer factorization up to one hundred decimal digits. For integers bigger then one hundred decimal digits  the general number field sieve is faster. We have not chosen to study the quadratic sieve because it seemed that it would take too much time to implement.



\section{Pollard-Rho Implementation}
We started by implementing Pollard's original algorithm \textsc{Pollard-Rho}(n) (algorithm \ref{Pollard-Rho}) in Java. That is the algorithm that uses Floyd's cycle detection.

After this we attempted to use Brent's cycle detection, which is an improved version of the Pollard's original algorithm called \textsc{Pollard-Rho-Brent}(n) (algorithm \ref{Pollard-Rho-Brent}). The algorithm was also implemented in Java. Unfortunately we were unable to improve the result, probably as an result of worse filtering of which numbers to factor.

Instead of trying to improve the previous implementation further we decided to quickly move on to the improvement of Pollard's algorithm by Brent, \textsc{Pollard-Brent}(n) (algorithm \ref{Pollard-Brent}), which is the best version of the algorithm that we could find. However we also wanted to try if we could improve our result by using a language other than Java. Due to the fact that none of us have had any programming experience in C++ (and only very limited in C) we decided to split the final programming effort into two parallel parts: one of us (Klaus Nicosia) started to implement \textsc{Pollard-Brent}(n) and the other two (David Falk and Marcus Öberg) started from scratch to implement \textsc{Pollard-Rho}(n) and then later \textsc{Pollard-Brent}(n) in C++.

Due to the time limits on the project we were unable to gain any significant improvement with our C++ implementations.

We achieved the best performance with our implementation of \textsc{Pollard-Brent}(n)  in Java for which the code is presented in appendix \ref{appendix_a} and appendix \ref{appendix_b}. The filtering strategies presented in section \ref{filtering} below concerns this implementation.

For all algorithms we used the function $f(x) = x^2 -1\;(mod\;n)$.

In Java we used \texttt{BigInteger} and in particular the methods \texttt{probablePrime} and \texttt{gcd} was used. Other common math functions from \texttt{BigInteger} and \texttt{Math} were used. No factoring method provided by Java (if any such exists) were used.

\subsubsection{Filtering the Numbers to Factorize}
\label{filtering}
We filtered what numbers to factor in two ways: by time limit and by number of iterations for different loops.

The \texttt{while}-loop at lines 90 through 106 in the method \texttt{pollardBrent} in \texttt{IntegerFactoringPollardRho.java} in appendix \ref{appendix_b} is where the program spends most of its time. Because of this the current time is recorded in this loop and compared to a time limit that is defined to be 14.5 s. If this time limit is exceeded the algorithm is aborted for the current number and \texttt{fail} is printed for this number and the following numbers that has not yet been factorized. This filter tries to maximize the amount time, of the allowed 15 s, that is used to factorize the one hundred numbers. Due to the random nature of the algorithm and the iteration limits described below this time limit filter is seldom enforced.

The most important iteration limit filter is also implemented in the \texttt{while}-loop in the \texttt{pollardBrent} method. This is what is called the \texttt{brentLimit} in the code and the loop breaks if \texttt{k} exceeds this number, note that \texttt{k} is increased by \texttt{m} at each iteration, were \texttt{m} is the number of bits of \texttt{n}. Through testing, a suitable value for \texttt{brentLimit} was found to be \texttt{32800} (or less for small \texttt{n}) for this problem. That is, the value \texttt{32800} is tailored specifically for the factoring of the one hundred numbers on the \textsc{Kattis} hardware.

Should the first part of the \textsc{Pollard-Brent}(n) method fail to calculate a factor of \texttt{n}, then \texttt{n} itself will be found. Thus the code in the \texttt{if}-statement at line 109 through 118 will run Pollard's original algorithm starting at the saved \texttt{ys} value that was stored before the batch \textsc{GCD} was calculated. The iteration limit for the loop in this part is \texttt{pollardLimit} and is at most \texttt{100}. Through testing we found that this loop seems to only iterate a few times (often less than ten iterations), therefore \texttt{pollardLimit} is chosen this small. We also found that there were seldom more than ten factors (of all the one hundred integers) for which Pollard's original algorithm had to be run, thus this iteration limit filter is almost never enforced.

As mentioned above  \texttt{m} is the number of bits of \texttt{n}. By testing different values for \texttt{m} we found around \texttt{100} worked best for our problem, that is to factor the integers which are close to one hundred bits or slightly lower.


\section{Results}
With \textsc{Pollard-Rho}(n) in java we were able to fully factorize 48 (\textsc{Kattis} ID 251926) out of the one hundred numbers and with C++ 49 out of one hundred numbers (\textsc{Kattis} ID 256012).

The \textsc{Pollard-Brent}(n) implementation in Java, combined with better filtering of which numbers to factor, was able to fully factorize 65 (\textsc{Kattis} ID 257621) out of the one hundred numbers to factor.
\section{Discussion}
When implementing \textsc{Pollard-Rho} in C++ we achieved a small improvement from our previous java implementation but we think that this improvement could have been larger if we had more time or more experience with C++. \textsc{Pollard-Brent} in java gave us a big improvement. 

When we profiled the code we saw that most of the CPU time was used by \texttt{GCD}. One of the reasons for the improvement with \textsc{Pollard-Brent} is probably because \textsc{Pollard-Brent} gets away with calling \textsc{GCD} much fewer times. We hoped to achieve an even better result by writing \textsc{Pollard-Brent} in C++ but this was not successful. Our C++ implementation got a runtime error in \textsc{Kattis}, a problem we were never able to recreate in out own test cases. Our guess was that the best choise of algorithm would have been the quadratic sieve, but we thought the algorithm was to complicated to implement in this short amount of time.

In out own tests we found that the different Pollard implementations had problems with numbers that was a power of a prime. This was solved by forcing the algorithm to calculate the x and y values through the f function at least one time with an do-while loop instead of an ordinary while, and to let it try again with a new random value if it failed. In general the hardest numbers to solve are those that consists of only two big prime factors. This was tested and found to work best with our C++ implementation of pollard-brent, but as previously mentioned this code was not accepted by \textsc{Kattis}.
\section{Conclusion} 
The Pollard-Brent algorithm was proved to be faster then the Pollard-Rho and Pollard-Rho-Brent algorithms.

\begin{thebibliography}{9}

\bibitem{cormen}
  Cormen and Leiserson and Rivest and Stein,
  Introduction to Algorithms,
  MIT Press,
  2nd Edition,
  2001.
  
\bibitem{wikipedia_pollard}
  Wikimedia Foundation (2011),
  Wikipedia, Pollard's rho algorithm,
  [www] <\url{http://en.wikipedia.org/wiki/Pollard's_rho}> 2011-11-16.

\bibitem{brent}
  Brent (1980),
  An Improved Monte Carlo Factorization Algorithm,
  [www] <\url{http://wwwmaths.anu.edu.au/~brent/pd/rpb051i.pdf}> 2011-11-16.

\end{thebibliography}

\appendix
\section{IntegerFactoring.java}
\label{appendix_a}
\lstinputlisting{../Java/src/IntegerFactoring.java}
\section{IntegerFactoringPollardRho.java}
\label{appendix_b}
\lstinputlisting{../Java/src/IntegerFactoringPollardRho.java}

\end{document}
